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THE  PSEUDO-SHOCK :  A  NON-LINE, Hi 
PROBLEM  OF  TRANSLATIONAL  RELAXATION 


B.  L.  Hicks 


Abstract 

The  Boltzmann  equation  hus  been  solved  by  Nordsieck' s  Monte  Carlo  method 
for  the  case  of  translational  relaxation  of  u  gas  of  elastic  spheres  whose  in¬ 
itial  velocity  distribution  function  has  the  form 

f(v,0)  =  ^|exp[ -n(v-lu)2]  +  exp[ -n (v+iu)2] j  . 

The  Mach  number 

M  =  (6n/5)^u 

describes  the  relative  separation  of  the  two  peaks  of  the  bimodal  distribution 
function  and  therefore  controls  the  degree  of  initial  departure  from  equilibri¬ 
um.  Calculations  have  been  made  for  M  =  0.5,  1,  2,  4,  ur.d  6,  which  includes  a 
range  of  initial  conditions  from  very  close  to  very  far  from  thermal 
equilibrium. 

In  the  early  part  of  the  translational  relaxation  we  find  that  relaxation 
by  a  factor  of  e  1  requires,  on  the  uveruge,  1.27  +  0.044  collisions  for  the 
lateral  temperature  and  0.80  +  O.O35  collisions  for  the  Boltzmann  function  H. 
The  temperature  relaxation  rote  is  thus  smaller  by  a  factor  of  I.58  +  0.120 
than  the  entropy  relaxation  rate.  (The  uncertainties  are  stated  as  90$  confi¬ 
dence  limits.)  These  collision  numbers  are  essentially  independent  of  M  and 
time,  it  least  until  each  molecule  has  made  about  ten  collisions.  Our  calcu¬ 
lations  agree  with  earlier,  more  qualitative  results  in  the  literature  that 
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correspond  to  different  initial  conditions.  We  have  also  shown  that  in  a 
Krook  model  of  our  relaxation  process,  the  ratio  of  the  two  collision  num¬ 
bers  is  somewhat  smaller  than  two  late  in  the  relaxation  and,  as  shown  al¬ 
so  by  Offerhaus,  approaches  two  asymptotically. 

Solution  of  the  Boltzmann  equation  for  any  other  reasonable  initial  con¬ 
ditions  could  be  obtained  with  the  same  Monte  Carlo  program.  For  example,  the 
relaxation  of  the  asymmetric  bi-modal  distribution  functions  used  by  Mott- 
Smith  to  describe  shock  wave  structure  could  be  studied. 
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1.  Introduction 

The  series  of  internal  CSL  reports,  "Numerical  Studies  of  Strong  Shock 
Waves",  '  describes  successive  steps  in  our  efforts  to  calculate  the  vel¬ 
ocity  distribution  functions  within  a  shock  wave  by  solving  the  Boltzmann 
equation.  We  have  shown  that  Nordsieck' s  Monte  Carlo  method  of  evaluating  the 
collision  integral  in  the  Boltzmann  equation  appears  to  be  both  reliable  and 
feasible.  Earlier  methods  of  evaluation  have  been  less  accurate  and  effici- 
ent.  3  The  chief  difficulty  yet  remaining  to  be  overcome  in  the  shock 

problem  is  that  of  insuring  an  adequate  rate  of  convergence  of  the  iterative 
method  of  solution  of  the  Boltzmann  equation.  We  therefore  decided  to  apply 
the  Monte  Curio  method  to  a  simpler  problem  than  the  shock  wave  in  order  to 


1.  Numerical  Studies  of  Strong  Shock  Waves.  Part  I.  Illiac  Solutions  of  a 

Boltzmann  Difference  Equation  by  Nordsieck' s  Method.  B.  L.  Hicks  and  J.  K. 
Aggurwul.  CSL  Report  1-111.  (1962) 

2.  _ .  Part  II.  Results  of  Illiac  Calculations.  B.  L.  Hicks  and 

J.  K.  Aggarwal.  CSL  Report  1-117.  (1963) 

3.  _ .  Part  III.  Studies  of  the  Monte  Carlo  und  Integration  Pro¬ 
grams.  B.  L.  Hicks.  CSL  Report  1-122.  (1963) 

4.  _ .  Ptrt  IV.  Description  of  the  1604  Program.  J.  K.  Aggarwal 

lid  B.  L.  Hicks.  CSL  Report  1-12}.  (1963) 

5.  .  Part  V.  Equations  und  Sculing.  B.  L.  Hicks.  CSL  Report 

I-12V.  (1963) 

6.  _ .  Part  VI.  Subroutines  and  Tables.  B.  L.  Hicks  und  M.  A. 

Smith.  CSL  Report  1-125*  (1964) 

7.  "Molecular  Dynamics  by  Electronic  Computers".  B.  J.  Alder  and  T.  Wain- 
wright  in  Transport  Processes  in  Statistical  Mechanics  (Interscience  Publish¬ 
ers,  1958),  PP*  97-131. 

8.  "Investigation  of  the  Many  Body-Problem  by  Electronic  Computers",  B.  J. 
Aider  and  T.  Wainwright,  in  The  Many  Body-Problem  (Interscience  Publishers, 
1963),  PP*  511-522. 


9.  G.  A.  Bird,  "Approach  to  Translational  Equilibrium  in  n  Rigid  Sphere  Gas". 
Phys.  of  FI.  6,  1518  (1963). 


demonstrate  at  once  the  utility  of  the  method  in  a  real  problem  and  to  learn 
from  this  exercise  some  new  tricks  to  use  on  the  shock  structure  problem. 
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We  chose  at  Nordsieck' s  suggestion  the  simpler  problem  of  calculating  the 
translational  relaxation  of  a  uniform  gas.  In  particular,  we  chose  as  the  in¬ 
itial  out -of- equilibrium  condition  a  bimodal  velocity  distribution.  This  con¬ 
dition  resembles  that  which  obtains  within  a  shock  wave,  and  could  be  made 
more  like  a  shock  wave  by  simple  extensions  of  the  calculations  reported  here. 
Hence  we  call  our  problem  of  relaxation  the  "pseudo-shock"  problem. 

We  have  solved  the  Boltzmann  equation  numerically  for  initial  conditions 
ranging  from  near  thermal  equilibrium  to  very  far  from  thermal  equilibrium. 

The  numerical  results  are  readily  interpretable  physically  and  suggest  a  new 
numerical  technique  that  should  be  useful  in  the  shock  wave  problem. 

We  are  indebted  to  Mrs.  Margaret  Smith  who  has  performed  most  of  the  pro¬ 
gramming  and  data  reduction  for  our  study  of  the  pseudo-shock. 

2.  Formulation  of  the  problem 

We  shall  treat  the  behavior  of  a  uniform  monatomic  gas  that  is  not  in 
thermal  equilibrium.  We  shall  consider  hard  sphere  molecules  of  diameter  a 
and  mass  m.  Molecules  with  other  force  fields  can  be  studied  if  their  differ¬ 
ential  cross  sections  are  known.  We  shall  suppose  that  there  are  no  external 
forces.  The  behavior  to  be  considered  is  then  that  of  a  simple  translationul 
reluxation  which  is  governed  by  the  Boltzmann  equation. 

One  set  of  units  will  be  used  in  the  body  of  the  paper.  A  second  set, 
"machine  units",  will  be  used  in  the  Sect.  4.3.  (See  also  Part  I(l).)  In 
the  first  set,  the  units  are  the  values,  denoted  by  the  subscript  1,  of  va¬ 
rious  properties  of  a  reference  gas.  Thus  n1;  T2  aie  the  units  of  number  den¬ 
sity  and  temperature.  The  unit  of  length  l1  =  l/(2nn1o2)  = 

((mean  free  path)  -y / -J2) .  The  unit  of  velocity 
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ci  =  (2nkT1/m)^  =  (mean  speed ) iX( n/2 ) .  The  unit  of  time  is  therefore 
(mean  free  time)!  and  of  the  velocity  distribution  function 
is  nj/cj3.  In  these  units  the  Boltzmann  equation  may  be  written 


df/dx  =  a  -  bf 


(FT*  -ff'  )|N*vrldv*  (dw/Un)  (2-1) 


where  f  =  f (v,x  )  is  the  time -dependent  velocity  distrib”tion  function,  x  is 
the  time  variable,  n  gives  the  direction  of  the  line  of  centers  during  a  col- 
lision,  v^  =  v  -  v,  and  f,  f  ,  F,  F  denote  the  four  values  of  f  correspond- 
ing  to  the  four  velocities  v,  v  ,  V,  V  characterizing  a  binary  collision. 

We  specify  as  the  initial  condition  at  X  =  0,  a  symmetrical  bimodal  vel¬ 
ocity  distribution  function  for  a  spatially  uniform  gas 

f(v,0)  =  fQ  =  ^jexp  [ -n(v-iu)2]  +  exp  [-n(v+iu)2]j  .  (2-2) 


At  all  later  times  the  velocity  distribution  function  depends  only  upon  x  and 
the  velocity  components  (v  ,v  )  which  are  cylindrical  polar  coordinates  in 
velocity  space.  In  our  units  the  parameter  u  is  related  to  a  Mach  number  by 
the  equation 

M  =  (6«/5)*u  .  (2-5) 

For  M  <  <  1  the  gas  described  by  fQ  is  almost  in  thermal  equilibrium  at  tem¬ 
perature  t  =  1  and  with  number  density  n  =  1.  For  M  >  >  1,  the  gas  described 
by  fQ  is  very  far  from  equilibrium,  and,  in  fact,  consists  of  two  oppositely 
directed  streams  of  gas,  each  at  temperature  t  =  1  moving  essentially  parallel 
to  the  x  axis  at  a  3peed  large  compared  to  the  thermal  velocity. 

We  shall  be  interested  in  the  behavior  of  the  velocity  distribution  func¬ 
tion  f,  the  collision  integral  (a-bf),  the  "lateral"  temperature 


and  the  Boltzmann  function 


H 


J 


f  log  f  dv  . 


(2-5) 


The  number  density  n  =  1  throughout  the  relaxation  process.  The  reduced 
temperature 

r  p 

T/Tj.  =  t  =  (2n/3n)  j  v2  f  dv=  l  +  ( 5/9 )Ma  =  1  +  ^  jiu2 

(2-6) 

is  likewise  constant  throughout  the  relaxation  process.* 

The  initial  and  final  (t  — >  00  )  values  of  tx  are  equal,  respectively,  to 


1  and  t.  The  initial  values 

of  M: 

of  H  can 

be  computed  analytically 

for  two  values 

H  ( 0 ) 

=  .  1 

2  * 

M  = 

0 

(2-7) 

H(0) 

=  .1  . 

2 

In  2  ,  M  = 

QD  . 

(2-8) 

For  other  values  of  M,  H(0)  must  be  computed  by  numerical  quadrature.  The 
asymptotic  value  of  H  is  giver,  by 


H(co  )  =  -  |  (1  +  In  t) 


(2-9) 


corresponding  to  the  asymptotic  velocity  distribution  function 

*00  (v)  =  =  t'  ^  exp  (-ttv2 /t )  . 


(2-10) 


*As  we  shall  see  later  (and  compensate  for)  n  and  t  do  not  remain  exactly 
constant  during  the  numerically  calculated  relaxation  process. 
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3.  Theoretical  Expectations 

From  Boltzmann' s  Theorem  we  know  that 

dfi/dx  <0  (3-1) 

during  the  relaxation.  The  asymptotic  value  of  H  for  an  equilibrium  gas  with 
given  t  and  for  n  =  1  is  given  by  Eq.  (2-9)>  and  according  to  Eq.  (3-1)  this 
value  must  be  the  lowest  value  reached  by  H  during  the  relaxation. 

Let  5f  =  6f(v,x )  represent  the  departure  of  the  velocity  distribution 
function  from  f  =  f(v,  oo),  its  equilibrium  value.  Then 

f  =  f  +  5f  (3-2) 

oo 

and 

r 

=  t  +  (fl/n)  j\’A2  6f  dv  .  (3-3) 

An  equation  similar  to  (3-3)  holds  for  any  moment  of  f. 

It  will  be  convenient  to  interpret  the  time  scale  and  collision  rates  in 
terms  of  the  behavior  of  the  equilibrium  gas,  described  by  f  .  We  may  expect 
the  time  constant  for  the  equilibrium  gas  and  for  the  asymptotic  part  of  the 
relaxation  to  be  inversely  proportional  to  the  mean  molecular  speed  and  there¬ 
fore  to  t  ^  .  The  asymptotic  relaxation  for  different  Mach  numbers  should 
then  be  similar  when  plotted  against  (t  a  t).  The  number  of  collisions  7V  suf¬ 
fered  by  a  molecule  in  the  equilibrium  reference  gas  in  (reduced)  time  Ax  is 

7L  =  (V2/n)t^  AX  =  0.450 1"5  At  .  (3-4) 

The  relaxation  will  be  almost  complete  when  (t^O  >  >  10.  In  this 
asymptotic  range,  j&f/f^l  ((  1  and 
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H 


H 


<■  ^  /  fif2  C  <!  -  Sf/3f  )  + 


OD 


ro 


00 


(3-5) 


^he  term  in  (6f)  is  missing  from  Eq.  (3~5)  because  n  and  t  are  constant  during 
the  relaxation.  Offerhaus^*^  has  discussed  the  (6f)2  term  but  not  the  (5f)3 
term  in  this  equation. 

To  get  a  more  detailed  view  of  the  relaxation  process  let  us  describe  the 
process  approximately  by  the  Krook  model,  in  which  the  Boltzmann  equation 
is  replaced  by 

df/ax  =  b(^-f).  (3-6) 


In  this  equation  If  is  a  Maxwell -Boltzmann  distribution  function  having  the 
same  (constant)  values  of  n  and  t  as  does  the  (non -equilibrium)  function  f. 
Therefore,  =  f^  is  independent  of  t.  The  solution  of  Eq.  (3-6)  then  is 


f 


f 

oo 


+ 


f  -f 
0  00 


(3-7) 


where  b  may  or  may  not  depend  upon  v,  according  to  the  nature  of  the  Krook 
model.  Notice  that  f,  in  Eq.  (3-7),  is  a  linear  combination  of  f  and  f 
We  may  thus  say  that  Eq.  (3-7)  corresponds,  for  our  translational  reluxation 
problem,  to  a  Mott-Smith  model  of  the  process  as  well  as  to  a  Krook  model. 

We  now  have  an  explicit  form  for  bf  to  use, 

6f  -  (f0-f00)e’bT  ;  (t*T  >>  1)  •  (3-8) 


10.  M.  F.  Offerhaus,  Theory  of  Relaxation  Phenomena  in  a  Monatomic  Gas.  Akad. 
64,  115  (1961),  see  Sect.  5 • 

11.  Bhatnaghar,  Gross,  and  Krook.  "A  Model  for  Collision  Processes  in  Gases." 
Fhys.  Rev.  a4,  511  (1954). 
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The  relaxation  of  tx  (or  other  moments)  is  similar  to  that  of  f: 


tx  =  t  +  (n/n)  J  (f0~f00)\2e~]:fZd^  • 


(3-9) 


But  for  large  enough  times,  namely,  for 


6f/3f  I  -  (f  -f  )e“DT/3f  «  1  , 

/>aol  OCD  '  ^  00  ^  1 


(3-lCj 


we  find  that 


K 


e 


-2oT 


dv 


(3-H) 


The  time  constant  b(v)  of  the  (asymptotic)  relaxation  for  each  velocity 
bin  increases  with  Increase  of  the  speed  v.  Let  b-L  be  the  smallest  value  of 
b.  Then  for  large  enough  times  the  relaxation  of  the  corresponding  velocity 
bin  will  dominate  the  relaxation  process  and  we  may  write  in  place  of  Eqs. 

(3-9),  (3-H) 

tj_  =  t  +  (l-t)e"biT  (3-12) 


II  =  H  +4  e"2blX  f  (f  -f  )2  f'1 


QD 


J 


O  OO  ' 


QD 


dv 


(3-13) 


Results  of  the  same  form  follow  from  the  alternative  assumption  that 
b  *  constant.  For  example,  we  may  assume  that  b2  in  the  Krook  model  corres¬ 
ponds  to  the  value  b(0)  for  the  equilibrium  gas;  ,hen 


bx  -  b(0)  -  t  */n  . 


(3-iM 


Insofar  as  the  Krook  model  is  valid  for  the  asymptotic  part  of  the  relax¬ 
ation  process  we  thus  can  predict  that  the  t ime  constant  b^  for  the  relaxation 
of  H  is  one-half  as  large  as  the  time  constant  b  =  b^  for  the  asymptotic 
relaxation  of  t  .  It  may  also  be  shown  for  the  Krook  model  that  the  term 
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in  Sf3  is  negative,  for  large  enough  values  of  M,  which  has  the  apparent  ef¬ 
fect  of  making  b^/bj.  ^  l/p_  when  b^T  ro  l.  Even  without  appealing  to  the  Krook 
model  we  can  show  that  this  is  a  reasonable  result  when  we  look  at  the 
(Bf)3/^  term  in  the  integrand  of  Eq.  (3-4).  This  term  behaves  near  a  resid- 

-  3  r 

ual  peak  like  t3  and  elsewhere  like  t  ~  .  Remembering  that  Bfdv  =  0,  we 
see  that  for  large  enough  t,  the  integral  of  Bf3/:^  is  positive  and  the  cor¬ 
responding  contribution  to  H  is  negative. 

For  large  values  of  M  we  can  also  make  a  prediction  about  the  early  part 
of  the  relaxation  process.  In  this  case  few  collisions  occur  except  between 
molecules  moving  with  velocities  u  and  (-u)  nearly  parallel  to  the  x  axis. 

For  t  ^  t  <<"  1  only  a  small  number  of  molecules  will  have  collided  and  these 
will  be  distributed  uniformly  in  velocity  space  along  the  circle  v  =  u.  We 
may  expect  then  (for  M  ^  )  1  and  (t  ^  T )  1)  that  we  will  see  in  isoline 

plots  of  f  a  "bridging"  along  a  circular  arc  between  the  two  delta  functions 

at  v  =  +u ,  v.  =  0. 

x  -  f  -b 

(12) 

It  should  be  noted  that  one  analytical  treatment'  '  of  translational  re¬ 
laxation  has  been  based  upon  superposition  throughout  uhe  relaxation  of  two 
Maxwe 11 -Boltzmann  distributions  the  separation  of  whose  centers  continuously 
decreases.  Th±  ■  model  cannot  account  for  the  directional  randomization  just 
discussed,  but  addition  of  a  third  Maxwell-Bolt zmunn  distribution  might  yield 
asymptotic  results  similar  to  those  suggested  by  the  Krook  model. 

We  are  now  ready  to  consider  the  numerical  methods  used  in  solving  the 
translational  relaxation  problem. 


12.  K.  Suchy,  "New  Methods  in  the  Kinetic  Theory  of  Rarified  Gases."  Ergeb. 
exakt.  Natural ss.  10}  (1964). 
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4.  Numerical  Methods 

4.1  Evaluation  of  the  Collision  Integral  by  Monte  Carlo  Sampling 

The  most  difficult  part  of  the  translational  relaxation  problem  is  evalu¬ 
ation  of  the  collision  integral.  It  is  fortunate,  sir^e  no  other  suitable 
method  (analytical  or  numerical)  is  yet  available,  that  Nordsieck's  Monte  Car¬ 
lo  method  is  well  developed. 

The  method  was  described  in  the  first  report  of  the  series.^7  The  col¬ 
lision  integral  is  first  replaced  by  an  integral  over  a  finite  region  of  vel¬ 
ocity  space  ohat  is  of  volume  R  and  that  includes  most  of  the  molecules.  The 
average  value  of  the  integrand  over  this  region  is  then  approximated  by  the 
average  of  a  large  and  fair  sample  of  N  values  of  the  integrand.  The  vulue  of 
the  collision  integral  is  then  given  by  the  product  of  this  average  value  with 
the  volume  R.  The  integrand  is  a  function  of  eight  independent  variables  de- 

-2*  -A 

rived  from  v;  v  ,  and  n.  Nordsieck1 s  Monte  Carlo  method  insures  fairness  of 
the  sampling  in  the  eight-dimensional  space  of  these  variables.  The  sampling 

and  Monte  Carlo  quadrature  make  use  of  226  velocity  bins  in  the  (v  ,v  )  plane. 

X  X 

In  a  numerical  solution  of  the  Boltzmann  equation,  it  is  the  speed  and 
accuracy  of  the  evaluation  of  the  collision  integral  that  must  be  our  primary 
concern  before  we  look  at  ether  characteristics  of  the  over-all  method  of 
solution.  For  samples  of  moderate  size  and  velocity  distributions  that  cover 
about  200  velocity  cells,  statistical  fluctuations  contribute  the  only  signi¬ 
ficant  error,  for  the  strtlsticc  1  error  is  then  much  larger  than  the  quadra¬ 
ture  error.  (Truncation  error  is  generally  small  except  for  large  values  of 
the  speed  v  where  the  collision  integral  itself  is  small.)  Thus,  with  our 
present  Monte  Carlo  program,*  the  calculation  of  226  values  of  (a-bf),  for  a 


*Detailed  tests  of  this  program  have  not  yet  been  described  in  reporti 
Tests  of  earlier  programs  are  described  in  Part  III. (5) 
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sample  of  N  =  104  collisions,  is  performed  in  50  sec.  on  the  CDC  1604  computer 
and  yields  statistical  errors  (expressed  as  probable  errors)  in  a  and  bf 
individually  of  +15 The  bias,  owing  to  quadrature  error,  amounts  to  +2 ‘jo  in 
uncorrected  a  or  bf  for  a  well-covered  velocity  space  for  an  equilibrium  gas. 
The  bias  in  (a-bf)  is  not  now  directly  determinable  for  a  non -equilibrium  gas 
but  is  reduced  by  the  method  described  in  Sec.  4.4.  We  would  like  to  empha¬ 
size  that  some  bias  in  (a-bf),  caused  by  quadrature  error,  is  to  be  expected 
whether  the  numerical  integration  uses  sampling  techniques  or  not. 

The  computing  time  is  proportional  to  N,  and  the  statistical  error  is 
proportional  to  II  5  .  For  many  calculations  it  is  not  practical  nor  is  it 
necessary  to  reduce  the  statistical  error  to  the  same  level  as  the  bias.  We 
have  obtained  significant  results  in  the  relaxation  problem  with  relatively 
small  samples  (N  -  104). 

For  M  ))  1,  only  a  few  velocity  bins,  in  effect,  are  used  in  Monte  Carlo 
(or  other)  numerical  evaluation  of  the  collision  integral.  The  large  result¬ 
ing  quadrature  error  then  presents  the  primary  difficulty  in  solving  the  re¬ 
laxation  problem  adequately  for  large  Mach  numbers.  (See  Sect.  6.32.) 


4.2  Integration  of  the  Boltzmann  Equation 


Consideration  of  the  various  sources  of  error  and  of  factors  influencing 
the  computing  speed  indicated  that  Euler  quadrature  in  the  time  variable  would  give  ap¬ 
propriate  accuracy  and  speed.  For  each  velocity  bjn  the  integration  formula  is 


f(t1+Ac)  =  f(Tx)  +  AT (df/dr ) 


(4-1) 


In  view  of  the  discussion  in  Section  3  we  chose  AT  to  be  approximately  propor¬ 
tional  to  t~  ^  so  that  quadrature  errors  in  the  forward  integration  are  ap¬ 
proximately  equal,  on  a  fractional  basis,  for  different  Mach  numbers. 


We  may  add  two  further  notes  about  the  forward  integration  process.  A 
process  of  higher  order  than  the  Euler  would  require  recalculation  of  (a-bf) 
one  or  more  times  before  a  step  in  time  is  made.  Improvements  in  accuracy  of 
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the  time-wise  integration  must  take  into  account,  however,  that  by  far  the  most 
time-consuming  part  of  the  whole  problem,  even  with  the  help  of  a  Monte  Carlo 
process,  is  the  evaluation  of  the  collision  integral.  Thus,  (see  Sect.  4.1) 
one  forward  step  in  time  tukes  50  sec.,  for  a  sample  of  104  collisions.  We 
can,  however,  on  the  basis  of  a  reasonable  assumption,  easily  correct  the  re¬ 
sults  of  the  Euler  integration  and  thereby  achieve  ubout  the  same  accuracy  in 
calculation  of  t  and  H  as  though  we  had  used  a  higher  order  integration 
process. 

The  assumption  used  is  that  each  function  (f,  t  and  H)  which  we  later 
discuss  is  nearly  an  exponential  function,  exp  (-pt^t).  For  the  function  f(t ), 
for  example,  we  may  then  write  the  expression  for  the  derivative  at  the  mid¬ 
point  of  an  interval  as 

f'  (T  -  -bt^  f(T  +  ^At)  £  -Pt^  f(T)  [1  -  ipAr]  . 

(4-2 ) 

Each  increment  Af  from  un  Euler  integration  should  then  be  multiplied  by  the 
(constant)  correction  factor  [l  -  ^f3t^At]  and  logarithmic  slopes  should  be  cor¬ 
rected  by  the  same  factor.  We  assume  that  the  logarithmic  slopes  of  the  func¬ 
tions  t^  and  H  should  be  corrected  by  the  factors  [l  -  if^t^AT]  and 

[l  -  ^•S./traAf]  where  f  and  are  determined  later.  The  correction  made 
2  H  x  H 

of  the  slopes  of  the  H  curves  is  greater  than  that  of  the  t  curves  but  is  less 
than  15$.  Higher  order  corrections  would  not  change  the  values  of  the  slopes 
by  more  than  1*. 
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4.3  Integrals  involving  f  and  (a-bf) 

( 1  5 ) 

As  in  our  earlier  work.  we  calculate  n,  nt,  nt  and  H  (eich  of  which 

involves  integration  of  f  over  velocity  space),  md  dn/dt ,  d(nt)/di  (each  of 
which  involves  integration  of  (a-bf)  over  velocity  space)  by  numerical  quad¬ 
ratures  in  velocity  space.  Combined  quadrature  and  truncation  errors  have 
been  reduced  to  less  than  1  ‘jo  for  all  values  of  M  for  large  values  of  (t  ) 
but,  for  large  Much  numbers,  (M  ^  4),  it  is  not  possible,  for  a  given  number 
of  velocity  bins  (226  in  our  c.se),  to  Keen  the  quadrature  errors  this  small 
for  ( t^T  )  <)  .  ( .rue  (  i  scission  f  Tale  TI  in  beet.  6.2.) 

To  f  cl.it  tc  control  of  the  trunc  tion  error  vc  introduced  Lwj  re!,  ted 
parameters  Kx  and  p.  The  par-meter  Kx  is  used  to  scale  the  variable  v  so  .s 
to  "fill"  the  velocity  space,  for  a  given  function  f.  Thus 


v 

m 


KjV 


(4-3) 


where  226  fixed  values  of  (vx;vL)[n  (velocity  components  in  "machine"  units) 
describe  the  velocity  cells  in  the  truncated  region  R  over  which  numerical 
integrations  are  carried  out.  Though  it  is  net  in  principle  necessary,  we 
keep  Kx  fixed  during  the  time-wise  integration  of  the  Boltzmann  equation  for 
any  one  value  of  M. 

As  a  simple  measure  of  truncation  error  (in  calculating  integrals  over  f 
we  use  the  integer  p  in  the  equation 

roD(2k/K1)  -  10-4  ;  (4-4) 

v  -  24  being  the  radius  of  the  spherical  region  R.  Combination  of  Eqs.  (4-4) 
and  (2-10)  then  gives  a  value  of  Ka  for  each  value  of  M  or  t .  We  took  p  to 
be  4 . 


i 
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We  have  not  examined  the  fractional  truncation  error  in  the  evaluation  of 
(a-bf).  We  would  expect  it  to  be  less  than  the  corresponding  errors  in  integ¬ 
rals  over  f  because  the  integrand  of  the  collision  integral  contains  products 
of  velocity  distribution  functions  and  therefore  decreases  much  more  rapidly 
than  f  us  the  speed  v  increases. 


4.4  Least  Square  adjustment  of  Calculated  Collision  Integrals 


As  v/e  have  noted  above,  some  bias  in  values  of  the  collision  integral 
calculated  by  numerical  quadrature  is  unavoidable.  To  reduce  the  resultant 
accumulation  of  errors  in  the  forward  integration  of  the  Boltzmann  equation  in 
time  we  have  devised  a  least  square  adjustment  of  the  collision  integrals  that 
are  calculated  by  Monte  Carlo  sampling.  This  method  is  preferable  to  the  "mo- 
ment  correction"  method'  ;  1  1  used  for  the  shock  problem  two  years  ago  be¬ 
cause  it  less  often  produces  negative  and  therefore  unacceptable  values  of  f. 

Let  us  first  consider  the  vulues  of  n  and  t  calculated  from  the  f 1  s  gen¬ 
erated  by  our  numerical  treatment  of  the  Boltzmann  equation.  These  values  of 
n  and  t  would  remain  constant  throughout  the  relaxation  process  if  there  were 
no  error  in  calculating  them.  We  will  go  far  toward  enforcing  constancy  of 
these  values  if  we  adjust  the  values  of  (a-bf)  as  little  as  possible,  in  a 
leust  square  sense,  while  imposing  the  two  conditions 


dn/dT  =  j  (df/dt  )dv  =  0 


(4-5) 


d(nt )/ dr 


(2*/3 


v^af/dT  )dv 


0 


(4-6 ) 


Both  the  Monte  Carlo  calculation  of  the  function  (df/dt)  and  the  numerical 
quadrature  used  to  approximate  the  integrals  in  Eq.  (4-5,6)  introduce  errors 
in  the  values  of  dn/dT  and  d(nt)/dT.  As  noted  before,  the  adjustment  process 
therefore  maintains  the  constancy  of  the  computed  values  of  n  und  of  t. 
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The  least  square  condition,  subject  to  Eqs.  (4-5)  and  (4-6),  is 

6  1  <VPs)2fsa  =  0  (4'7) 

^0 

in  which  p  and  P  are  the  values  of  df/dt  (for  the  sth  velocity  bin)  before 
s  s 

and  after  adjustment.  Solution  of  the  least  square  problem  yields  a  simple 

explicit  formula  for  P  in  terms  of  moments  of  df/dt  that  are  approximated  by 

s 

weighted  sums.  Note  that  we  minimize  the  mean  square  value  of  the  f ractlonal 
adjustment  of  p  ,  u  procedure  that  is  consistent  with  a  characteristic  of 
Nordsieck'  s  Monte  Carlo  method,  namely,  the  approximately  constant  frictional 
error  of  the  values  of  df/dt  that  it  produces. 

We  checked  the  effectiveness  of  the  correction  method  by  obtaining  Monte 
Carlo  solutions  of  the  Boltzmann  equation  for  M  =  0  both  with  and  without  the 
correction  of  (n-bf).  For  M  -  0  the  gus  is  initially  in  equilibrium  so  that 
properties  of  the  gas  that  are  calculated  by  the  Monte  Carlo  method  should  de¬ 
viate  from  the  equilibrium  values  of  these  jroperties  only  because  of  fluctua¬ 
tions  ;nd  bias  in  the  Monte  Carlo  method. 

With  a  program  which  uses  the  correction  method,  the  f  and  a  isolines  for 
M  =  0  differ  in  only  a  few  bins  from  the  equilibrium  isolines.  This  is  true 
both  after  two  and  20  steps  in  time.  The  isoai.ues  of  (a-bf)  indicate  probable 
erro  ~s  of  about  24$relative  to  a.  When  no  corrections  of  (a-bf)  have  been 
made  then  t,  t  and  |h|  increase  by  1.6,  1.6  and  1.7 1°,  respectively,  after  the 
20  steps  in  time.  Correction  of  (a-bf)  yields  of  course  a  constant  computed 
value  of  t  (low  by  0 . 3%  owing  to  quadrature  und  truncation  error  in  integrat¬ 
ing  f  to  get  t)  md  yields  values  of  t  and  H  that  fluctuate  near  the  equili¬ 


brium  values.  These  fluctuations  in  amount  to  ~0.4^£>  and  in  H  to 'vO. 05/6. 
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The  mean  value  of  |Hl  as  calculated  from  the  corrected  (a-bf)  is  low  by  O.yjo 
(quadrature  and  truncation  error  over  f)  and  perhaps  shows  a  slight  downward 
bias  amounting  to  0. 1%  in  20  steps  in  time.  The  size  of  the  fluctuations  will 
be  used  later  in  interpreting  the  pseudo-shock  calculations. 

These  results  show  that  the  least  square  adjustment  of  the  values  of 
(a-bf)  does  indeed  reduce  the  trends  uway  from  equilibrium  that  are  produced 
by  unavoidable  slight  bias  in  the  Monte  Carlo  (or  other)  calculation  of 
(a-bf). 

4.q  Monte  Carlo  Fluctuations 

As  in  an  earlier  investigation^  we  wished  to  study  statistical  varia¬ 
tion  of  quantities  that  had  been  calculated  with  the  help  of  Monte  Carlo  sam¬ 
pling.  Therefore  we  made  four  runs  of  the  same  type  but  each  run  with  a 
different  and  independent  sample.  For  each  run  we  chose  to  use  a  sample  of 
N  -  104  collisions  for  each  of  ten  steps  in  time,  starting  at  T  =  0  for  each 
Mach  number,  and  also  at  T  =  1CAT ,  so  that  we  could  judge  the  effects  of  the 
Mnch  number  and  of  the  phase  of  the  relaxation  upon  the  variances  among  the 
runs . 

We  shall  be  concerned  chiefly  in  later  sections  with  the  variance  of  the 
dcriv  tives  of  t  and  H  as  functions  of  (t^t). 


5.  Pa  rarnerers 

The  parameters  that  define  the  numerical  treatment  of  the  translational 
relaxation  process  are  At,  N,  p,  and  the  first  random  number  used  in  generat¬ 
ing  sequence  of  independent  samples  for  a  given  run.  Except  us  specified 
otherwise,  the  sequences  of  samples  were  wholly  independent  of  one  another. 
There  is  only  one  physical  parameter  governing  the  relaxation  process,  n  r.e- 
ly  M.  F  rameters  used  in  the  various  runs  re  summarized  in  Table  I. 


M 
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Table  I.  Summary  of  Pseudo  Shock  Calculations 


M 

interval 

AC 

t^  AT 

Range  of 

(t/ac) 

£ 

no. 

runs 

of 

N 

0 

.25 

.25 

1-20 

4 

1 

10,923 

1-20 

1  (b) 

•  5 

.25 

.2663 

1-10 

4 

4 

10,923 

11-16 

1 

11-30 

1 

l 

.25 

.5119 

1-10 

4 

4 

10,923 

11-20 

4 

11-30 

1 

2 

.25 

.4488 

1-10 

4 

4 

10,923 

11-20 

4 

11-30 

1 

.125 

.2244 

1-10 

4 

1 

43,691 

4 

.125 

0931 

1-10 

4 

4 

10,923 

11-20 

4 

11-30 

1 

.0625 

.1965 

1-4 

4 

1 

43,691 

6 

.0625 

.2864 

1-20 

4 

1  (c 

) 

10,923 

1-20 

3 

1  (c 

) 

1-20 

2 

1  (c 

) 

1-10 

4 

3 

11-20 

3 

21-30 

1 

10 

.0625 

.4700 

1-10 

b 

4 

10,923 

11-20 

4 

21-30 

1 

(a)  All  runs  were  made  with  independent  samples  unless  indicated  otherwise. 

(b)  No  (u-bf)  corrections. 

(c)  These  three  runs  were  made  with  the  same  sample. 
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6.  Results  and  Discussion 


Preliminai 


irison  with  Theoretical 


;ctations 


In  interpreting  our  computational  results  we  shall  be  interested  in  the 
behavior  of  the  lateral  temperature  tj_  and  the  Boltzmann  function  H,  and,  to 
a  lesser  extent,  in  the  detailed  appearance  of  isoline  plots  of  f,  a,  and 
(a-bf).  The  behavior  of  t^  and  H  is  shown  in  Figs.  1  and  2  for  five  values  of 
the  Mach  number  M.  In  each  figure,  the  abscissa  is  t  where  t  is  the  (re¬ 
duced)  equilibrium  temperature  given  by  Eq.  (2-6)  and  t  is  the  (reduced)  time 
variable.*  The  ordinates  in  Figs.  1  and  2  represent  the  absolute  value  of  the 
differences  (t^(0)-t^M(t ) )  and  (H' (oo)  -  H^t),  calculated  from  the 
Monte  Carlo  results,  where  t^(0)  and  H' (®  )  are  the  asymptotic  values  of  tj_ 
and  H.  Discussion  of  the  estimation  of  these  and  other  asymptotic  values  will 
be  deferred  to  later  sections.  Let  us  now  state  briefly  what  these  figures 
show  that  was  predicted  in  Sect.  3* 

The  numerical  results  do  satisfy  Boltzmann's  Theorem,  at  least  until  the 
relaxing  gas  is  close  enough  to  equilibrium  so  that  the  Monte  Carlo  fluctua¬ 
tions  produce  fluctuations  of  H  above  and  below  the  equilibrium  value  H'(oo). 
Except  for  M  -  0.?,  the  r  nge  of  (H.j(f)  -  no))  before  such  deviations  occur  is 
two  or  more  decades.  Both  tj_  and  H  (which  are  measures  of  departure  froi:. 
equilibrium)  relax  as  though  a  single  relaxation  mechanism  were  effective  from 
the  onset  of  the  relaxation  process.  Time  constants  for  different  values  of  M 

are  proportional  to  +  ^  ,  s  indicated  by  the  parallelism  of  the  relaxation 

X 

curves  plotted  with  t  T  on  the  abscissa.  The  time  constant  for  relaxation 


*We  recall  that  0. 


is  the  number  of  collisions  that  have  occurred 


in  reduced  time  T  in  the  reference  gas.  (See  Eq.  (3-4).) 


of  H  does  appear  to  be  ubout  one-half  as  large  as  for  the  relaxation  of 
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In  Fig.  3  the  predicted  "bridging"  is  clearly  visible  in  the  f-isolines 
for  M  =  4.  These  isolines  (see  Table  I)  were  calculated  for  AT  =  0.0625  cor¬ 
responding  to  u  collision  of  only  one  molecule  in  11  in  the  reference  gas, 
showing  the  sensitivity  of  Nordsieck1  s  method.  Our  examination  of  the  iso¬ 
lines  for  f,  a  and  (a-bf)  also  shows  that  they  appro. ich,  for  t  T  ))  1,  the 
shapes  characteristic  of  equilibrium,  as  we  would  expect. 

Although  a  more  careful  discussion  of  these  points  will  be  given  later, 
we  can  already  see  that  agreement  with  the  theoretical  expectations  gives  con 
siderable  indirect  support  to  the  validity  cf  our  Monte  Carlo  solution  of  the 
Boltzmann  equation  for  the  pseudo-shock. 


6.2  Asymptotic  Behavior 

Let  us  first  examine  the  behavior  of  tj_  and  H  for  large  times,  i.e.  for 
T  =  20  AT.  We  have  made  a  number  of  different  estimates,  given  in  Table  II, 
of  both  the  initial  and  asymptotic  values  of  t  and  H.  These  estimates  are 
derived  from  different  combinations  of  analytical  formulas,  numerical  integ¬ 
rations,  and  Monte  Carlo  calculations.  It  is  necessary  to  consider  various 
estimates  of  n(t )  and  t(t)  as  well  as  of  t±  (t  )  and  H(c). 

Let  us  use  the  subscript  a  to  indicate  a  wholly  analytical  calculation; 

the  subscript  q  to  indicate  result  derived  by  numerical  integration  of  m 

u  _ 

analytical  f(v, c)  over  v;  and  the  subscript  M  to  indicate  a  result  derived  by 
numerical  integration  of  a  Monte  Curio  f (v,T )  over  v.  Initial  and  asymptotic 
values  of  the  various  estimates  are  indicated  by  the  values  0  and  oo  of  the 


argument  r  . 
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Figure  J.  Early  stages  in  the  relaxation  of  the  velocity  distribution  func¬ 
tion  f  for  a  Mach  number  of  4.0.  v^  and  v  are  cylindric  .1  polai 

coordinates  in  velocity  space.  The  (reduced)  time  interv  il 
At  -  0.0625  corresponds  to  a  collision  of  one  molecule  in  eleven 
in  a  reference  gas. 
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Table  II.  Limiting  Values  of  Macroscopic  Properties  of  the  Relaxing  Gas 


1. 

2. 

3- 

4. 

5- 

6 . 

7- 

8. 

O. 

i0. 

M 

n 

u 

n  (0) 

q 

n  (00  ) 

Q 

t 

a 

t  (0) 
q 

V°°> 

t  (0) 
j.q 

t  (00 ) 

jV  ' 

w°° 

0.0 

1.000 

■99973 

1.000 

.99743 

.99736 

0.5 

1.000 

. 99982 

.99956 

1.1339 

1.1363 

1.1334 

.99317 

1.1333 

1.14 

1.0 

1.000 

1.0002 

•99995 

1.5556 

1.5539 

1.5500 

.99893 

1.5499 

1.55 

2.0 

1.000 

1.0015 

1.0012 

3.2222 

3.2201 

3.2119 

.99691 

3.2117 

3.22 

4.0 

1.000 

1.0148 

1.0146 

9 . 3389 

9.3696 

9.8436 

.96952 

9.3429 

9.87 

6 .  o 

1.000 

1 . 0730 

1.0777 

21.000 

20.965 

20.313 

.83621 

20.812 

20.3 

10.0 

1.000 

1.2914 

1.2913 

56.556 

53.976 

53.876 

.92132 

53.871 

oo 

1.000 

00 

1. 

11. 

12. 

13. 

14. 

15. 

16  . 

M 

H  (o) 

a 

H  (0) 

q 

H  (00 ) 

a 

H  (oo  ) 

4. 

H1  (00  ) 

v®  > 

0.0 

-1.5000 

-1.4953 

-1.5000 

- 1 . 4  96 1 

0.5 

- 1.6696 

-1.6951 

-1.6373 

-1.6915 

-1.639 

1.0 

-1.9557 

-2.1623 

-2.1570 

-2.1614 

-2.160 

2.0 

-2.1797 

-3-2551 

-5.252-3 

-3-2575 

-3.254 

4.0 

-2.1947 

-4.9371 

-4.9971 

-4.9923 

-4.989 

6.0 

-2.1361 

-6.0663 

-6.4430 

-6.4436 

-6.442 

.0.0 

-2.8734 

-7.5523 

-9-3284 

-9.3331 

oo 

-2.1532 

-00 

vn 
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In  the  second  column  of  Table  II  appear  the  analytical  values  of  n(x), 
namely  n^,  which  has  the  value  unity  at  all  times  because  the  number  of  mol¬ 
ecules  is  conserved.  In  column  3  is  the  estimate  n^(0),  obtained  by  numerical 
integration  over  the  analytical  values  of  f(v,0)  given  by  Eq.  (2-2).  Trunca¬ 
tion  errors  are  negative  for  n  (0)  (as  they  are  also  for  numerical  Integra- 
tions  yielding  estimates  of  j  H(x  )  |  and  of  all  the  positive  functions  in 
Table  II)  and  decrease  as  M  increases.  The  quadrature  errors  in  n^(0)  are  ap¬ 
parently  positive  and  increase  (unavoidably,  for  a  fixed  number  of  velocity 
bins)  as  M  increases.  The  error  in  t^(0)  (column  6)  is  negative  and  smaller 
in  absolute  value  than  the  error  in  n^(0).  (Note  that  numerical  quadrature 
yields  values  of  (nt)  and  (nt^)  rather  than  t  and  t^  directly.) 

We  must  emphasize  at  this  point  that  our  digital  computer  solution  of  the 
Boltzmann  equation  produces  and  uses  values  of  tj_,  H,  dn/dT  etc.  that  have  been 
obtained  by  numerical  integration.  In  "clamping"  the  values  of  n(T )  and  t(x), 
by  the  method  described  in  Section  4.4,  we  are  then  fixing  the  values  of 
VT>’  t^(x  )  at  their  initial  values,  n^(0),  t^(0).  As  the  translational  re¬ 
laxation  proceeds,  the  quadrature  errors  (over  v)  in  H  and  t^  decrease  and  are 
much  less  than  for  all  Mach  numbers  for  (t  )  ))  1.  (The  quadrature  er¬ 

rors  in  (a-bf )  presumably  decrease  even  more  rapidly..  )  The  proper  interpreta¬ 
tion  of  our  relaxation  calculations  is  then  that  we  are  following  the 


relaxation  of  a  gas  whose 


sotic  density  and  temperature  are  n^(0)  and 


t^(0) .  We  shall  therefore  base  several  estimates  of  the  asymptotic  values  of 
other  functions  upon  these  values  n^(0)  and  t^(0). 

The  velocity  distribution  function  for  a  gas  in  thermal  equilibrium  and 
having  the  density  and  temperature  n^(0),  t^(0)  is 


fQ0(v)  =  n^(0)[ t^(0)] ”  *  exp  [ -flv^/t^O)]  . 


(6-1) 
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(Compare  eq.  (2-10).)  Numerical  integration  over  this  function  gives  the 
values  n^(ac),  t^(co  )  listed  in  columns  4  and  7  of  Table  II.  The  differ¬ 
ences  between  n^(0)  and  n^(oo  )  and  between  t^(0)  and  t^(oo)  are  very  small, 
which  testifies  to  the  small  combined  truncation  and  quadrature  errors  of  our 
numerical  integration  over  v  for  functions  f(v)  that  "fill"  properly  the  part 
of  v  space  we  use  for  numerical  integration. 

Having  given  this  much  detailed  discussion  of  the  quantities  appearing  in 
the  first  seven  columns  of  Table  II,  we  can  discuss  rather  briefly  the  related 
quantities  listed  in  columns  8,  9>  11,  12,  13,  14,  15  of  the  Table.  The  ini¬ 
tial  values  of  t  (0)  *n  column  8  exhibit  small  deviations  from  the  correct 

value  t  *1  for  low  Mach  numbers  and  larger  errors  for  high  Mach  numbers 

la 

much  as  did  t^(0).  Also,  as  we  might  expect,  t  (oo  )  in  column  9  agrees  very 
well  with  t  (oo)  for  all  values  of  M. 

q 

There  seems  to  be  no  analytical  formula  for  H  (0)  as  a  function  of  M. 

The  initial  and  asymptotic  (M  — *  oo )  values  shown  in  the  Table  in  column  11 
are,  however,  given  by  the  Eqs.  (2-7>8)*  Note  the  small  range  of  H  as  com¬ 
pared  to  that  of  t.  The  initial  values  H^(0)  obtained  by  numerical  Integra- 
tion  over  v  are  given  in  column  12  and  show  large  departures  from  the  correct 
values  only  for  M  =  10.  Exact  values  of  H  ( oo )  (column  13)  are  given  by  Eq. 

3. 

(2-9)  and  may  be  compared  with  the  values  of  H^(oo)  (column  14)  obtained  by 
numerical  integration  of  [f(v,oo)ln  f(v,oo)]  over  v.  The  agreement  is  good  ex¬ 
cept,  as  for  M  =  6,  10,  where  n  (0)  is  much  different  from  unity.  In  view  of 

Q. 

our  earlier  discussion  a  more  useful  comparison  is  that  between  H^(oo)  and 


H'(oo)  =  n  (0)[ln  n  (0)  -  1.5  In  t^(0)  -  1.5] 


(6-2) 


As  we  would  expect,  these  two  sets  of  values  of  H(oo  )  agree  within  much  less 
than  1$. 
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We  are  now  ready  to  discuss  the  behavior  of  the  rel  :xing  gas  s  deter¬ 
mined  from  the  Monte  Carlo  calculations.  From  our  discussion  in  Section  J 
we  expect  that  t  (t)  and  K(t }  will  each,  to  first  approximation,  relax  ex¬ 
ponentially  to  their  equilibrium  values.  To  explore  this  possibility  we 

1  , 

plotted  [t  (0)  -  t^^(T)j  and  lH^(r)  -  H'(oo)]  on  logarithmic  scales  vs.  (t  T 
i  linear  scule  as  in  Figures  1  md  2.  For  analysis  of  the  early  part  of  the 
relaxation  process  the  asymptotic  values  chosen  for  subtraction  from  t.^(c 
and  H^(T )  need  be  vulues  only  approximately  equal  to  t  (0)  and  H'(oo).  For 
the  later  parts  of  the  relaxation,  where  the  departure  from  the  equilibrium 
values  is  small,  it  is  necessary  to  choose  these  values  to  be  t  (0)  and  H' (oo  ) 
(or  other  values  such  as  t  (oo)  and  H^(oo)  that  approximate  closely  to  them). 

The  straight  und  parallel  parts  of  the  curves  in  Figures  1  and  2  will  be 

discussed  in  u  later  section;  we  are  concerned  here  only  with  the  appurent 

asymptotic  behavior*  of  the  functions  t^(t  )  and  From  these  cur/es  we 

derived  estimates  of  the  asymptotic  values  of  it  (0)  -  t^^(oo)]  und  of 

[H^(oo)  -  H' (oo  )]  and  from  these  calculated  values  of  t^  ( oo )  and  H^(oo  )  as 

gi'ven  in  columns  10  and  1 6  of  Table  II.  The  values  of  t  (oo)  .re  equal, 

within  the  'uncertainty  of  esti.. rating  them  (say  O.yjo),  to  the  values  of  t  (0). 

The  values  of  H^(oo)  for  M  4  10,  are  equal  to  the  values  of  H^(oo  )  or  H'(oc) 

within  0.2/o  or  less,  an  error  level  comparable  with  the  difference  between 

H  (oo  )  und  H1  (oo  ) . 

0. 


^Careful  numerical  analysis  of  the  Krook  model  shows  that  we  would  not  be 
able,  in  the  Monte  Carlo  calculations,  to  find  at  what  time  the  quantity 
[Hj^(c)  -  H'(oo)]  becomes  accurate?  r  proportional  to  it  (0)  -  t^(  r  )] 2,  say 

within  2/b,  because  of  the  smallness  of  these  differences  and  the  very  gradual 
change  of  slope  of  the  curves  in  the  asymptotic  region. 


on 
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We  nay  thus  conclude  that  oar  Monte  Carlo  integration  of  the  Boltzmann 
equation  produces  the  correct  usynptotic  behavior  of  t^  and  H  except  for  small 
errors  that  are  to  be  expected  from  the  nature  of  the  numerical  integrations 
over  v. 

We  may  make  a  few  remarks  about  the  asymptotic  statistical  fluctuations 
of  tj^(t )  and  of  H^(r).  The  fluctuations  of  t^,  (for  M  y  0)  are  roughly  equal 

to  the  fluctuations  in  t  found  in  Section  4.4  for  an  initiul  state  of  equi- 

J.M 

librium  (M  =  0).  The  fluctuations  of  H^(r)  (for  M  0),  however,  amount  to 
about  0.01^  for  all  Mach  numbers  and  -are  therefore  smaller  by  a  factor  of  about  five 
than  were  found  for  M  -  0.  Note  that  this  remarkably  small  level  of  fluctuation  was 
obtained  with  a  rather  small  Monte  Carlo  sample,  namely,  104  collisions. 

6.3-  Behavior  for  (t  ^  10 

6.31  m  4  6 

Two  striking  features  of  the  carves  in  Figures  1  and  2  ure  evident  for 
M  ^  6:  they  are  almost  straight  and  are  nearly  parallel.  To  discover  how 
nearly  straight  and  parallel  we  can  suy  they  are,  we  need  to  examine  in  detail 
various  sources  of  error.  We  restrict  our  first  discussion  to  the  first  ten 
steps  in  At  (for  which  (t  3)* 

Data  to  permit  estimation  of  the  fluctuations  in  slope  owing  to  the  Monte 
Carlo  calculations  of  the  collision  Integral  were  obtained  by  the  method  de¬ 
scribed  in  Section  4.5*  A  typical  "fan"  generated  by  this  method  is  shown  in 
Figure  4.  Since  we  were  primarily  interested  in  the  apparent  constancy  of  the 
slope  us  M  changes  we  decided  to  examine  first  the  mean  slope  and  the  (sample) 
standard  deviation  S  for  each  set  of  four  runs  of  ten  intervals  each.  The 
slopes  were  determined  by  fitting  a  straight  line  (that  passed  through  the 


4 


Figure  4. 


i 


_ I _ L. 

10  1  5 

Collision  number 


Temperature  relaxation  computed  with  four  independent  seqcenc 
Monte  Carlo  samples.  Each  sample  contains  10,92}  collisions, 
number  M  -  2.0. 
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point  for  x  =  0)  to  the  curve  for  each  sample.  The  results  are  given  in  Ta¬ 
ble  III  and  hove  been  corrected  for  quadrature  error  in  integrations  with  re¬ 
spect  to  X.  (See  Section  4.1.)  It  is  unnecessary  to  correct  for  the 
difference  between  tQ  and  t  (0)  in  the  abscissa  variable  (t^x )  because  this 
difference  is  so  small.  (See  Table  II.) 

Inspection  shows  that  the  variation  among  the  mean  slopes  for  the  differ¬ 
ent  Mach  numbers  is  generally  no  larger  than  the  90^  confidence  limits  for 
each  Mach  number.  To  make  a  more  careful  judgement  we  take  the  following 
steps : 

( 13 ) 

a)  Comparison  of  S  -  from  variance  ratio  tests'  J '  we  conclude  that  the 
sample  standard  deviations  for  the  deferent  Mach  numbers  are  not  significant¬ 
ly  different  at  the  5$  level. 

b)  Comparison  of  the  means  -  an  analysis  of  variations^)  shows  that 
the  sample  means  for  the  slopes  of  the  K  curves  are  not  significantly  differ¬ 
ent  from  one  another  at  the  5$  level.  The  mean  slopes  for  t^  agree  with  one 
another  even  better  than  do  the  mean  slopes  for  H. 

We  may  thus  conclude  that  the  logarithmic  slopes  of  the  relaxation  curves 
(plotted  v£  (t  ) )  are  independent  of  M  for  t  5  and  M  ^  6,  and  are 

equal,  respectively,  to  0.3539  +  0.0122  and  O.56O7  +  O.O233,  the  mean  values 
of  the  slopes  for  t^  und  H.  These  values  of  mean  slope  and  confidence 
limits  agree  with  values  obtained  less  formally  for  the  range  5  4  ^ 

Detailed  analysis  in  this  range  is  less  profitable  because  the  estimated  slopes 
are  sensitive  (for  (t  )  near  10)  to  statistical  fluctuations,  to  the  exact 
choice  of  the  nominal  asymptotic  values  of  the  two  variables,  and  to  quadra- 
ture  and  truncation  errors  in  integratiuns  over  v. 


13.  u.  E.  Weatherbum,  A  First  Course  in  Mathematical  Statistics  (Cambridge 
lQ6l)  Chap.  X  and  XI. 
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Table  III. 

Monte  Carlo 

Fluctuations 

of  Slope  of 

Relaxation 

Curves 

*l 

H 

M 

mean 

slope 

sample  s.d. 

S 

90* 

confidence 

limits 

mean 

slope 

sample  s.d. 
S 

9C % 

confidence 

limits 

0.5 

.3598 

.0139 

.0189 

.5596 

.0338 

.0459 

1.0 

0451 

.0133 

.0180 

.6120 

.0197 

.0267 

2.0 

0536 

.0150 

.0203 

.5638 

.0268 

.0363 

2.0* 

.3528* 

.0150* 

. 0203* 

.5948* 

.0268* 

.  0363* 

4.0 

•  3599 

.0430 

.0585 

•  53^8 

.0842 

.1143 

6.0 

•  3522 

.0565 

•  0767 

.4989 

.0708 

.0961 

mean 

•  3539 

.5607 

*The  starred  values  of  mean  slopes  were  obtained  from  one  special  run 
with  a  sample  four  times  larger  than  any  sample  used  for  the  other  runs.  For 
this  special  run  we  take  the  value  of  S  and  of  the  confidence  limits  to  be  the 
same  as  for  the  values  given  in  the  line  above  for  the  same  Mach  number  M. 


r 
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The  analysis  of  the  early  part  of  the  relaxation  process  hus  thus  far 
dealt  with  the  apparent  paraxlelism  of  the  curves,  each  curve  having  been  as¬ 
sumed  to  be  straight.  We  now  ask,  "How  straight  are  the  curves?"  Inspection 
of  the  mean  curves  for  each  fan  (find  also  the  curves  for  the  special  large 
sample  at  M  =  2.0)  show  no  indication  of  curvature  that  cannot  be  attributed 
to  statistical  fluctuations  or  to  small  uncertainties  in  the  choice  of  the 
asymptotic  values  of  t  (0)  and  H'(oo). 

Our  over-all  conclusions  are  that  t .  and  H  (for  M  ^  6  and  (t^r  <  10) 
each  relaxes  as  though  it  were  governed  by  a  single  relaxation  process. 

t±(T  )  =  1  +  (t-1)  exp  ( -f3x t^r  )  (6-6) 

H(t)  =  H(0)  +  lH(oo)  -  H(0)]  exp  (-PRt^r)  (6-7) 

where 

P  =  0.5539  +  0.0122  (6-8) 

0H  =  0.5607  +  0.0233  (6-9) 

and  ~e  independent  of  M.  The  uncertainties  are  stated  in  terms  of  cjOj>  confi¬ 
dence  limits  and  are  bused  upon  observed  Monte  Carlo  fluctuations.  The  uncer¬ 
tainties  do  not  include  estimates  of  possible  systematic  errors. 

The  ratio  (3„/f3  -  1.534  +  0.120  (90^  confidence  limit)  which,  as  pre- 

dieted  in  Sect.  3 >  Is  somewhat  smaller  than  the  known  asymptotic  value  of 
two. 

Eq.  (3-I3)  for  the  Krook  model  (bx  -c.  b(0))  corresponds  to 
6,  =  n  1  =  O.3I8.  We  may  interpret  the  lurger  value,  0.354,  obtained 


by  the  Monte  Carlo  solution  of  the  Boltzmann  equation  by  finding  for  what  mo¬ 
lecular  speed  b(v)  (for  the  equilibrium  gas)  is  equal  to  0.354  t ^  .  This 
speed  is  O.^lv  where  v  is  the  rr.ear.  molecular  speed.  Molecules  whose  speed  is 
less  than  0.51v  relax  more  slowly  than  those  whose  speed  is  larger  than  O.Slv. 
This  conclusion  could  be  examined  in  detail  by  analysis  of  the  time  variation 
of  the  velocity  distribution  functions  output  in  our  calculations  for  each 
velocity  bin. 

The  number  of  collisions  in  a  reference  gas  necessary  for  relaxation  by  a 
factor  of  e  1  is  V2/(jt(3)  =  0.450/(3.  The  relaxation  of  t^  and  of  H  may  then 
be  described  as  follows:  1.27  +  0.044  collisions  of  each  molecule  are  re¬ 
quired  for  relaxation  of  t^  by  the  fuctor  e  1 }  and  0.00  +  0,033  collisions  are 

required  for  relaxation  of  H  by  the  same  fuctor.  (The  provisional  value  of 

(7 ) 

this  second  number  obtained  by  Alder  and  Wainwright'  for  a  different  start¬ 
ing  condition  was  about  one.)  Each  number  is  representative  of  the  transla¬ 
tional  relaxation  process  initially  and  up  to  the  time  that  relaxation  by  a 
factor  of  c  2  or  e  3  has  occurred.  Asyraptot ically;  only  0.64  collisions  per 
molecule  may  be  necessary  for  relaxation  of  H  by  the  factor  e  .  (See 
Sect .  3.  ) 

6.32  M  >  6 

For  large  enough  values  of  M  we  snould  expect  that  the  computational 
methods  used  here  would  break  down.  Errors  then  would  enter  because  of  the 
small  number  of  velocity  bins  in  which  f(v, 0)  is  sensibly  different  from  zero. 
The  resulting  quadrature  errors  amount  to  29^  and  5 1°,  for  n  and  t  respectively, 
for  M  -  10.  Even  more  serious  quadrature  errors  arise  in  the  Monte  Carlo  (or 
any  other)  numerical  evaluation  of  the  collision  integral.  In  spite  of  these 
problems,  t^(t )  — ^  t  (0)  and  H^(t  )  — >  H'(oo)  for  M  =  10  as  T  increases.  How¬ 
ever,  the  behavior  of  t^,(T  )  end  H^(  1  )  for  small  times  exhibits  very  large 
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fluctuations.  We  have  not,  therefore,  presented  our  results  in  detail  for 
M  =  10. 

It  is  possible  that  semi- analytical  calculations  like  those  of  Alder  and 

( 7 ) 

Wainwright  would  give  accurately  the  detailed  pattern  of  translational  re¬ 
laxation  for  M  )  )  1.  Our  methods  could  be  used  in  relaxation  calculations 
for  values  of  M)  6  if  ve  used  a  substantially  larger  number  of  velocity  bins. 


